function posts=initial_posts(ests,b_mu,w_mu)

options=struct('Display','off','MaxFunEvals',100000,'MaxIter',500,'GradObj','off','TolFun',1e-5);

L=length(ests)/2;
pi_hat=ests(1:L,:);
delta_hat=ests(L+1:2*L,:);

gamma=zeros(L,2);
kappa=zeros(L,2);
mu=[b_mu,w_mu];
for l=1:L
    for col=1:2
        psi_hat=fsolve(@(x)fit_psi(x,pi_hat(l,col),delta_hat(l,col),mu(1,col)),[0;0],options);
        if abs(psi_hat(1))<3
            gamma(l,col)=psi_hat(1);
        end
        if abs(psi_hat(2))<3
            kappa(l,col)=psi_hat(2);
        end
    end
end

posts=[gamma(:,1);kappa(:,1);gamma(:,2);kappa(:,2)];